﻿//Credit to Paul Bourke (pbourke@swin.edu.au) for the original Fortran 77 Program :))
//Converted to a standalone C# 2.0 library by Morten Nielsen (www.iter.dk)
//Check out: http://astronomy.swin.edu.au/~pbourke/terrain/triangulate/
//You can use this code however you like providing the above credits remain in tact

using System;
using System.Collections.Generic;

using SharpGL.Media3D;
using SharpGL.Triangulator.Geometry;

namespace SharpGL.Triangulator
{
    /// <summary>
    /// Performs the Delauney triangulation on a set of points.
    /// </summary>
    /// <remarks>
    /// Based on Paul Bourke's "An Algorithm for Interpolating Irregularly-Spaced Data
    /// with Applications in Terrain Modelling"
    /// http://astronomy.swin.edu.au/~pbourke/modelling/triangulate/
    /// </remarks>
    public class Delauney
    {
        #region Private Methods

        /// <summary>
        /// Returns a triangle that encompasses all triangulation points.
        /// </summary>
        /// <param name="points">A list of triangulation points.</param>
        /// <returns>Returns a triangle that encompasses all triangulation points.</returns>
        private static Triangle FindSuperTriangle(IList<Point3D> points)
        {
            if (points == null && points.Count >= 1)
                throw new ArgumentException("points");

            // Find the maximum and minimum vertex bounds.
            // This is to allow calculation of the bounding supertriangle
            float xmin = points[0].X;
            float ymin = points[0].Y;
            float xmax = xmin;
            float ymax = ymin;
            for (int i = 1; i < points.Count; i++)
            {
                if (points[i].X < xmin) xmin = points[i].X;
                else if (points[i].X > xmax) xmax = points[i].X;
                if (points[i].Y < ymin) ymin = points[i].Y;
                else if (points[i].Y > ymax) ymax = points[i].Y;
            }

            float dx = xmax - xmin;
            float dy = ymax - ymin;
            float dmax = (dx > dy) ? dx : dy;

            float xmid = (xmax + xmin) * 0.5f;
            float ymid = (ymax + ymin) * 0.5f;

            // Set up the supertriangle
            // This is a triangle which encompasses all the sample points.
            return new Triangle(new Point3D((xmid - 2 * dmax), (ymid - dmax), 0.0f),
                                new Point3D(xmid, (ymid + 2 * dmax), 0.0f),
                                new Point3D((xmid + 2 * dmax), (ymid - dmax), 0.0f));

            /*double M = points[0].X;

            // get the extremal x and y coordinates
            for (int i = 1; i < points.Count; i++)
            {
                double xAbs = Math.Abs(points[i].X);
                double yAbs = Math.Abs(points[i].Y);
                if (xAbs > M) M = xAbs;
                if (yAbs > M) M = yAbs;
            }

            // make a triangle
            Geometry.Point sp1 = new Geometry.Point(10 * M, 0, 0);
            Geometry.Point sp2 = new Geometry.Point(0, 10 * M, 0);
            Geometry.Point sp3 = new Geometry.Point(-10 * M, -10 * M, 0);

            return new Geometry.Triangle(sp1, sp2, sp3);*/
        }

        #endregion Private Methods

        #region Public Methods

        /// <summary>
        /// Performs Delauney triangulation on a set of points.
        /// </summary>
        /// <remarks>
        /// <para>
        /// The triangulation doesn't support multiple points with the same planar location.
        /// points-lists with duplicate points may result in strange triangulation with intersecting edges.
        /// To avoid adding multiple points to your vertex-list you can use the following anonymous predicate
        /// method:
        /// <code>
        /// if(!points.Exists(delegate(Triangulator.Geometry.Point point) { return pNew.Equals2D(point); }))
        ///		points.Add(pNew);
        /// </code>
        /// </para>
        /// <para>The triangulation algorithm may be described in pseudo-code as follows:
        /// <code>
        /// subroutine Triangulate
        /// input : vertex list
        /// output : triangle list
        ///    initialize the triangle list
        ///    determine the supertriangle
        ///    add supertriangle points to the end of the vertex list
        ///    add the supertriangle to the triangle list
        ///    for each sample point in the vertex list
        ///       initialize the edge buffer
        ///       for each triangle currently in the triangle list
        ///          calculate the triangle circumcircle center and radius
        ///          if the point lies in the triangle circumcircle then
        ///             add the three triangle edges to the edge buffer
        ///             remove the triangle from the triangle list
        ///          endif
        ///       endfor
        ///       delete all doubly specified edges from the edge buffer
        ///          this leaves the edges of the enclosing polygon only
        ///       add to the triangle list all triangles formed between the point 
        ///          and the edges of the enclosing polygon
        ///    endfor
        ///    remove any triangles from the triangle list that use the supertriangle points
        ///    remove the supertriangle points from the vertex list
        /// end
        /// </code>
        /// </para>
        /// </remarks>
        /// <param name="points">List of points to triangulate.</param>
        /// <returns>triangles referencing vertex indices arranged in clockwise order</returns>
        public static List<Triangle> Triangulate(IList<Point3D> points)
        {
            if (3 > points.Count)
                throw new ArgumentException("Need at least three points for triangulation.");

            int trimax = 4 * points.Count;

            List<Edge> edgesBuff = new List<Edge>();
            List<Triangle> triangles = new List<Triangle>();

            // The supertriangle coordinates are added to the end of the
            // vertex list. The supertriangle is the first triangle in
            // the triangle list.
            Triangle superTriangle = FindSuperTriangle(points);
            // superTriangle placed at index 0
            triangles.Add(superTriangle);

            // Include each point one at a time into the existing mesh
            for (int i = 0; i < points.Count; i++)
            {
                // Set up the edge buffer.
                // If the point (points[i].X,points[i].Y) lies inside the circumcircle then the
                // three edges of that triangle are added to the edge buffer and the triangle is removed from list.
                for (int j = 0; j < triangles.Count; j++)
                {
                    if (triangles[j].ContainsInCircumcircle(points[i]))
                    {
                        edgesBuff.Add(new Edge(triangles[j].FirstPoint, triangles[j].SecondPoint));
                        edgesBuff.Add(new Edge(triangles[j].SecondPoint, triangles[j].ThirdPoint));
                        edgesBuff.Add(new Edge(triangles[j].ThirdPoint, triangles[j].FirstPoint));
                        triangles.RemoveAt(j);
                        j--;
                    }
                }

                // Remove duplicate edges
                // Note: if all triangles are specified anticlockwise then all
                // interior edges are opposite pointing in direction.
                for (int j = edgesBuff.Count - 2; j >= 0; j--)
                {
                    for (int k = edgesBuff.Count - 1; k >= j + 1; k--)
                    {
                        if (edgesBuff[j].Equals(edgesBuff[k]))
                        {
                            edgesBuff.RemoveAt(k);
                            edgesBuff.RemoveAt(j);
                            k--;
                        }
                    }
                }
                // Form new triangles for the current point
                // Skipping over any tagged edges.
                // All edges are arranged in clockwise order.
                for (int j = 0; j < edgesBuff.Count; j++)
                {
                    if (triangles.Count >= trimax)
                        throw new ApplicationException("Exceeded maximum edges.");
                    triangles.Add(new Triangle(edgesBuff[j].FirstPoint, edgesBuff[j].SecondPoint, points[i]));
                }
                edgesBuff.Clear();
            }

            // Remove triangles with supertriangle points
            for (int i = triangles.Count - 1; i >= 0; i--)
            {
                if (triangles[i].SharesVertexWith(superTriangle))
                    triangles.RemoveAt(i);
            }

            triangles.TrimExcess();
            return triangles;
        }

        #endregion Public Methods
    }
}